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We have studied the exchange-correlation hole n xc and pair correlation function in the valence 
shell of the Si atom in its 3 P ground state, using highly accurate Slater- Jastrow wavefunctions and 
the Variational Monte Carlo method. The exchange-correlation hole shows a number of interesting 
features caused by the open shell structure of Si, including a marked transition from efficient to poor 
screening behavior as a test majority-spin electron is moved from the center of a valence p orbital 
onto the axis perpendicular to the occupied p orbitals. This behavior results from the dramatic 
difference in the exchange hole in the two cases, which is partially compensated by a corresponding 
anisotropy in the correlation hole. In addition we observe an anisotropic change in the spin density 
induced by Coulomb correlation, reducing the spatial overlap between the different spin-components 
of the density and contributing to the anisotropy of the correlation hole. The exclusion effect 
correlation, in which a 3s 2 — * 3p 2 substitution excitation is allowed along the unoccupied axis and 
prohibited along the other axes, was found to have a noticeable effect on the correlation hole and is 
partially accountable for its anisotropy; however, it is inconsistent with the observed changes in spin 
density. In contrast to the longer range features, we find that the "on-top" correlation hole is well 
described by linear density functional theory, for a large range of local density and magnetization. 

PACS numbers: 31.25.Eb, 31.10.-t-z, 71.15.-m, 02.70.Lq 



I. INTRODUCTION 

An understanding of the exchange-correlation hole and 
related quantities such as the pair correlation function 
and exchange-correlation energy is an important factor 
in the systematic development of accurate density func- 
tional theory (DFT) methods for quantum chemistry and 
solid state physics [Q. Much progress has been made in 
understanding the pair correlation function in the homo- 
geneous electron gas or jellium, using accurate numeri- 
cal techniques |^|| and analytic modeling Q] . Although 
the closely related Coulomb hole has long been studied 
in atomic and molecular systems [^]-[7j, the quantitative 
understanding of correlation holes in inhomogeneous sys- 
tems is far from complete. In recent years there have been 
several efforts to fill the gaps, with a focus on character- 
izing properties closely tied to the development of den- 
sity functional theories, such as the on-top hole |q M| a nd 
system-averaged exchange and correlation holes [||fl0| . 

In much of the previous work on the exchange- 
correlation hole in atoms and molecules, the focus has 
been on two electron systems and closed-shell atoms such 
as Be and Ne, with additional studies of more compli- 
cated molecules such as N2 and H2O j7j. One class of 
systems that has received less attention is that of open- 
shell atoms. The valence shell of such an atom may con- 
sist only of a few electrons about a radially symmetric 
core, but with a degenerate ground state and the cor- 
responding valence-shell structure, have interesting and 



quite complex features in its exchange-correlation hole. 
In particular, in the case of less-than-half filling, these 
include the exclusion-effect correlations involving three 
or more electrons [ pd| . Modeling such systems requires 
proper attention both to orbital correlations best treated 
in a configuration interaction (CI) context and to dy- 
namic correlations more amenable to a density functional 
description. The application of density functional the- 
ory to open-shell systems is an active area of research, 
with the question of the optimal treatment of degenerate 
ground states and of the formation of such atoms into 
molecules giving rise to intriguing problems Jl^ j . 

Recently, Variational Monte Carlo (VMC) meth- 
ods combining highly accurate trial wavefunctions with 
Monte Carlo methods for evaluating ground state ener- 
gies and wavefunctions have been developed and applied 
to atoms and molecules Jl^-[l5| ■ The wavefunctions used 
in these calculations have the advantage of being simple 
and compact, typically recovering 85% or more of the 
correlation energy of atoms and the dissociation ener- 
gies of molecules with a few trial parameters. Unlike 
configuration interaction expansions they can describe 
with equal ease both "nondynamic" correlations and "dy- 
namic" correlations such as the short-range cusp condi- 
tion jl6j. These features make this method a natural 
candidate for studying electron correlations, and it has 
been employed in recent studies of the pair correlation 
function in crystalline Si jT7| , |l8| ]. 

In this paper we study the exchange-correlation hole 
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and pair correlation function of the valence shell of the 
Si atom with the VMC method. A quantitative knowl- 
edge of correlations in Si is important in improving DFT 
predictions for the cohesive energy, binding energies and 
surface characteristics for this and related technologically 
important materials. The Si atom is in its own right a 
useful laboratory for the study of electron correlations, 
with both local jellium-like and the nonlocal chemical 
properties due to the open shell structure of the atom 
playing important roles in determining the relevant phys- 
ical and chemical properties. The paper is organized as 
follows: Sec. || provides theoretical background on the 
exchange-correlation hole, Sec. Ill gives details of the sys- 



tem studied and of the calculational method. Results for 
the exchange-correlation hole are given in Sec. IV , for 
the pair correlation function in Sec. \V\ a nd for correla- 
tion effects on the spin density in Sec. [VI We close with 



a summary of the results and conclusion in Sec. VII 



II. THE EXCHANGE CORRELATION HOLE 

The exchange-correlation hole, including explicit spin 
dependence, is defined as: 



n xc (ro,a ;T,a) 



n (2) (r ,(7o;r,cr) 
n(r Q ,a ) 



n(r,a). (1) 



Physically it describes the change in density at r and 
for spin component a from its mean value, n(r, a), given 
the presence of another electron with spin cr at position 
r . The quantity n( 2 )(r , cr Q ; r, a) is the pair density, the 
expectation in the ground state of finding a pair of elec- 
trons at two given coordinates. It is defined in terms of 
an expectation of the ground state as: 

n (2) (r ,cr ;r, a) = V] (S{r - r i )(5 (ToCTi (5(r - Yj)8 acTj ). 



(2) 



A closely related quantity, the pair correlation function 
<7(ro, (To, r i a )i is a measure of the pair density relative to 
that expected for uncorrclatcd electrons with the same 
density distribution: 



;(r ,CT ;r,cr) = 



n (2) (r ,a ;r,g) 
n(r , <r )ra(r, a) ' 



(3) 



The importance of the exchange-correlation hole to den- 
sity functional theory lies in its connection jl9],|o) to the 
exchange-correlation energy E xc : 



E *c = \J d 3 r n(r) J d 3 r' J dX 



i f 1 ^ n xc (r,r',X) 



(4) 



Here n xc (r,r',X) is the exchange-correlation hole, 
summed over spins, for the system with scaled Coulomb 
interaction Ae 2 , with an external potential altered so that 



the density of the system remains unchanged. The in- 
tegration over coupling constant strength A accounts for 
the kinetic energy cost of correlating electrons, weakening 
the strength of the integrated n xc with respect to its value 
at A = l. Although E xc is in principle determined from a 
knowledge of the single particle density alone, this depen- 
dence is in general not easy to determine beyond the local 
density approximation. The exchange-correlation hole 
has a wealth of features which may be used to test and 
improve theoretical models of the exchange-correlation 
energy. As a result, many attempts to systematically 
improve density functional theory have this function as 
a starting point p]pl]-p3|. 

In this paper, we will discuss the full coupling-constant 
(A = 1) case for n xc . Although the coupling-constant inte- 
grated quantity is most directly connected to the density 
functional theory, the full coupling-constant case is inter- 
esting in itself, as its average is an experimentally mea- 
surable expectation of the ground state 0|. It is also 
an essential ingredient in modern hybrid methods which 
combine elements of density functional theories and con- 
ventional Hartree-Fock methods [ p2[ . The quantity n xc 
is often analyzed by a decomposition into an exchange 
component, n x , corresponding to the A = or noninter- 
acting system, which describes the correlations between 
particles arising from the Pauli exclusion principle, and 
a correlation component, n c , determined by taking the 
difference between the fully interacting and noninteract- 
ing cases, which describes the additional correlation due 
to the Coulomb interaction between electrons j2j|. We 
study the explicit spin decomposition of n xc . This choice 
is useful for understanding the correlation response of an 
open-shell atom for which the ground state has nontrivial 
differences in its spin components. Spin decomposition 
falls roughly along the lines of exchange and correlation: 
exchange affects only particles with the same spin, and 
Coulomb correlations, though also present in the same- 
spin case, are most noticeable in the opposite-spin chan- 
nel. 

Some initial insight into the nature of the exchange- 
correlation hole in atoms and other finite systems can be 
obtained by considering limiting cases. The exchange- 
correlation hole about a reference electron in the ho- 
mogeneous electron gas is isotropic, that is, a function 
solely of the distance from the electron, and localized 
about the electron with a radius determined by the aver- 
age interelectron distance. Consequently, in a system of 
slowly varying density, the the local density approxima- 
tion (LDA) holds, in which the hole is determined by the 
density (or each spin-component of the density if these 
are different) at the location of the reference electron and 
"moves" with the position of the reference electron. In 
atoms or molecules, important correlation effects often 
involve a pair excitation from the noninteracting ground- 
state into a finite number of lowlying, perhaps nearly de- 
generate, excited states. In this case, the resulting corre- 
lation hole is dependent on the shape of the ground-state 
orbitals vacated and the excited-state orbitals occupied 
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and is largely unsensitive to the position of the reference 
electron. In real systems, the exchange-correlation hole 
will contain aspects of both limiting cases, with those 
that "move" with the position of the electron termed 
"dynamic" correlations and the orbital correlations in- 
sensitive to electron position termed "nondynamic" . In 
open-shell atoms both play important roles. 



III. CALCULATION APPROACH 

A. Correlated description of the Si atom 

We focus on correlations in the valence shell of 
the atom with the core electrons replaced by norm- 
conserving ab initio nonlocal pseudopotentials derived 
from LDA calculations [^6| . In this case the valence shell 
of the atom is described by the Hamiltonian 



2m 



Vi i (r l )\km i ><li 



li ,7Ui 



E lr= 



r i -r.i 



(5) 



where the sums are over the valence electrons. Vi is the 
nonlocal pseudopotential and \livrn >< Urriil the single 
particle projection operator onto the state with total an- 
gular momentum and z-axis projection m;. The final 
term is the intravalence Coulomb interaction. 

The Si atom has a nine- fold degenerate (3s 2 3p 2 ) 3 P 
ground state. By maximizing the spin projection, this 
state can be represented by a single Slater determinant, 
consisting of a majority-spin component, here chosen to 
be spin up, with one 3s and two 3p orbitals, and a mi- 
nority or down-spin component with one 3s electron. 
There then remain two angular momentum projections 
that lead to physically significant differences in n xc . The 
mi = projection has p x and p y orbitals and a "pancake" 
like shape, while the m; = ±l projections have a "cigar" 
shape, with pg and p± orbitals. This distinction does not 
play a role in determining E xc for an atom, given the in- 
variance of the energy to rotations of the atom. The pres- 
ence of a quantization axis in the formation of a Si bond 
breaks this invariance, so that the projection-specific be- 
havior of n xc becomes important in determining accurate 
molecular binding energies. It thus should be useful as a 
test of density functional theory, which typically overes- 
timates binding energies by about 1 eV [QJlfj] . Results in 
this paper focus on the m/ = projection which provides 
a clear comparison between the situation perpendicular 
to and parallel to the p orbitals, though calculations were 
done on the other projection for comparison. 

To calculate the exchange-correlation hole variation- 
ally we start with a Slater-Jastrow trial wavefunction 



iff = exp(-F)Y[D a 



with D a being a Slater determinant for spin component 
a and F a Jastrow correlation factor. All single-particle 
orbitals are obtained from the same local DFT program 
that determined our pseudopotentials. We use a Boys 
and Handy form mA 27 for F, which includes electron- 



electron, electron-nucleus and electron-electron-nucleus 
correlations expanded in a basis set of correlation func- 
tions: 



F 



E 



Cj r . 



i*3 



{r\r 
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-m—l 



f l )f n 



(7) 



bri/(l + bri) and r. 



(6) 



The basis functions are f. 
drij/(l + drij), where is the distance between a pair 
of electrons, is the distance between electron i and the 
atom center. The terms 6, d and Q mn are variational 
parameters. The lowest order fy term is set separately 
for opposite and same-spin electron correlations to satisfy 
the short range electron-electron cusp condition jl6| for 
each case. Higher order terms treat longer range effects 
and are determined without distinguishing electron spin. 
Electron-ion terms in the correlation function correct for 
the tendency of interelectron correlations to expand the 
volume of the atom and provide density-dependent cor- 
rections to the electron-electron correlation |14|. Since 
the valence shell of Si is less than half-filled, orbital or 
nondynamic correlations may be important These 
can be incorporated with a multideterminant extension 
of the Slater-Jastrow wavefunction: 

<P = exp(-F)J2v«l[D%' (8) 



B. Method of Calculation 

The Variational Monte Carlo method |||] is used to 
calculate the ground-state energy, derivatives with re- 
spect to variational parameters, and other expectation 
values. The heart of the method lies in the judicious 
statistical sampling of integration points to obtain an 
estimate of the many-body integrals involved in evaluat- 
ing expectations of the Slater-Jastrow trial wavefunction. 
This estimate is limited in accuracy by statistical noise; 
however if the trial wavefunction is a good approximation 
to the ground-state wavefunction this noise can be very 
easily managed with a relatively small number of con- 
figurations |13J. Optimized wavefunctions are obtained 
by minimizing the variance of the energy With a 

trial wavefunction of 18 expansion terms including two 
set by the same- and opposite- spin cusp conditions, up to 
l + m + n = 6 in the basis function expansion, we obtain a 
value of 3.8028(2) a.u. for the ground-state valence-shell 
energy. The correlation energy, measured with respect 
to a noninteracting ground-state energy of 3.7188 a.u., 
is 97% of the correlation energy obtained from Green's 
function Monte Carlo and 95% of that obtained from CI 
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using the same nonlocal pseudopotential [IL5|. A simi- 
lar calculation starting from a two determinant reference 
point [Eq. (§h], adding the 3pl3p x 3p v excited state to the 
nonintcracting ground-state configuration, resulted in a 
modest improvement in energy to 3.8041(2) a.u. or 96.6% 
of the valence-shell correlation energy with respect to CI. 

To calculate correlation functions, we measure spin- 
decomposed single-particle densities n(r, a) and condi- 
tional densities n(r, <r|ro, To). The spin-dependent con- 
ditional density is defined as the ground-state density 
distribution as a function of spin a and position r of 
the N — 1 other particles given one with spin <7o fixed 
at ro. The difference between the conditional and "un- 
restricted" densities, gives the spin-dependent exchange- 
correlation hole, Eq. ([[]), 



n xc {r , ct ; r, a) = n(r, a\r , a ) - n(r, a) 



(9) 



Separate calculations are done to measure the density 
and the conditional density for various values of ro and 
o"o- These expectations are first calculated exactly for 
the Slater-determinant wavefunction (setting F = in 
our trial wavefunction). Then, the difference between 
the expectations obtained with the Slater determinant 
and the fully interacting wavefunctions is measured sta- 
tistically using Monte Carlo sampling and the method 
of correlated estimates 28 . This technique is an effi- 



cient means to estimate statistically the change in the 
expectation of an observable under a small perturbation 
of the Hamiltonian or of the variational parameters of 
the wavefunction, taking advantage of the high degree 
of correlation between the two expectations to reduce 
noise in the difference of their statistical estimates. In 
the present case, the correlation hole, which describes the 
difference between the correlated and Slater-determinant 
exchange-correlation holes, is fairly small compared to 
the exchange-correlation hole. It is thus a reasonable as- 
sumption that the correlated estimation approach should 
improve the sampling efficiency for this quantity. In prac- 
tice, the procedure was observed to reduce statistical er- 
rors in our data by a factor of 5, and roughly 10 5 random 
samples of the wavefunction sufficed to obtain expecta- 
tion values. 

The expectations for density and conditional density 
were expanded in a plane wave basis, taking the aver- 
age of exp(—iG ■ Ti) for a set of plane waves up to a 
32 Ry cutoff, on a supercell 18 Bohr radii (as) in length. 
The coefficients were symmetrized and fast- Fourier trans- 
formed to obtain real-space densities. This expansion 
provides smooth profiles despite the statistical noise of 
the sampling. On the other hand the short-wavelength 
cutoff causes an unrealistic rounding off of the correla- 
tion hole cusp at short interparticle distances, and spu- 
rious oscillations at low densities. At present this cutoff 
is the largest error in our calculation. In the case of the 
exchange-correlation hole this error is notable mostly in 
the cusp region. A more complete discussion of errors 
is presented in Sec. [v| in regard to the pair correlation 



function which is more sensitive to the cutoff error than 
the exchange-correlation hole. 



IV. RESULTS FOR THE 
EXCHANGE-CORRELATION HOLE 

The exchange hole in DFT is obtained by evaluating 
n xc from Eq. (yj) for the Slater determinant wavefunction 
that minimizes the energy in the noninteracting (A = 0) 
limit: with the Coulomb interaction replaced by a single- 
particle potential that reproduces the true ground-state 
density [^0) . In practice the Slater determinant of local 
DFT orbitals used in our calculations produces a density 
that differs from our VMC density by a few percent. In 
terms of these orbitals, the exchange hole is: 



n x (r ,a ;r, a) 



n(r ,a ) 



(10) 



where the sum runs over all the occupied single-particle 
orbitals tp a for the spin component a. This expression, 
as a function of r for fixed ro (i.e., interpreting the hole 
as the change in density of the system at r given a parti- 
cle observed at r ) has the form of minus the probability 
density of a hybrid atomic orbital. That is, it describes a 
normalized linear combination of single-particle orbitals 
with coefficients c a = tpa( r o, &o) / \/ n i r o , °o)- This choice 
of Cq, for each orbital ipa represents the unique linear 
combination of orbitals that maximizes the probability 
for an electron to be observed at ro and spin ctq (con- 
versely giving zero likelihood for any other electron to be 
observed at that point.) Finally the integral over r of the 
exchange hole is —1, as it measures the integrated differ- 
ence in density between the (N— l)-electron system given 
one electron fixed at ro and the full A^-electron system. 

For a spin-up particle on the z axis in the L z = pro- 
jection, this exchange-hole orbital is a 3s state since the 
occupied p states (p x ,p y ) are in the x-y plane. The ex- 
change hole is completely insensitive to the electron po- 
sition on this axis, given only one possible orbital from 
which to construct it. The situation along this axis thus 
corresponds to an extreme departure from the "dynamic" 
picture of the exchange hole derived from the homoge- 
neous electron gas. For a particle on the x axis, the hole 
is a combination of s and p x orbitals. We find that the 
occupation probability |c p | 2 of the p orbital to be between 
60% and 75% in the region (0.5 clb<x<2.5 a B ) of peak 
density along the x axis; this is roughly equivalent to that 
of a sp 2 hybrid orbital (which is two-thirds p) oriented in 
the direction of the fixed particle. The slight variation in 
the c p acts to keep the hole more or less centered on the 
reference electron in the region of peak density on the 
x-y plane. The exchange hole as a result is more sensi- 
tive to the exact location of particles in this plane, and 
therefore more efficient in screening them, than along the 
z axis, leading to significant differences in the correlation 
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FIG. 1. The same spin (a), opposite spin (b) and total (c) exchange-correlation hole about a spin-up particle located at 1.4 as 
from the atom center on the x axis (parallel to the p orbitals) for the ground state of the Si atom in the L z — projection. The 
surface plot shows the change in density along a plane cutting through the origin along the x and z axes. 




holes in the two cases as well. This transition from "non- 
dynamic" screening along the axis perpendicular to the 
occupied p orbitals to a more "dynamic" screening in the 
plane occupied by them, particularly as it occurs at peak 
density in the valence shell, constitutes a major difference 
between the n xc of Si and that of closed-shell atoms. 

In Figs. [I] and we plot the exchange-correlation hole 
about a spin- up electron fixed at the point of peak density 
parallel to and perpendicular to the two 3p orbitals of the 
mi = (pancake) projection of the Si atom. The single- 
configuration Slater- Jastrow wavefunction [Eq. (^)] is 
used. Each plot shows the response to this electron in the 
x-z plane, that is, the plane cutting through the center 
of the atom at the plot origin, with one axis (x) parallel 
to and one (z) perpendicular to the occupied p orbitals. 
The hole is split into same spin (a), opposite spin (b) 
and total (c) response. The comparison between these 
two situations shows the dramatic anisotropy in n xc re- 
flecting that of the exchange hole. 

For a spin-up electron placed on the x axis, Fig. fi], the 
exchange-correlation hole is dominated by the sjr-like 
exchange hole. The additional effects of Coulomb corre- 
lation on the same-spin channel are hard to detect, while 
the opposite-spin correlation hole is small (contributing 



14% of the total on-top hole, or value of the hole at zero 
interparticle separation.) It is largely confined to a nar- 
row region about the electron, indicating that the sp 2 
hybrid hole screens the electron efficiently. 

Fig. ^ shows the exchange-correlation hole of a spin- 
up electron on the z axis, perpendicular to the two p 
orbitals. With the addition of correlation, the same-spin 
hole (a) loses the rotational symmetry of the 3s state that 
characterizes the exchange hole, with the polarization of 
the two Sp orbitals creating a double valley on either 
side of the reference electron. The opposite-spin hole (b) 
shows the polarization of the 3s spin-down orbital, with a 
well centered about the fixed electron and a strong dipole 
response at longer range. 

The total z-axis exchange-correlation hole, (c), shows a 
smooth interpolation of the two spin contributions lead- 
ing to a large crescent-shaped area near the electron 
from which the other electrons are repelled. The correla- 
tion hole contributes considerably to the total exchange- 
correlation hole, with up to 40% of the total on-top hole 
due to correlation. Although the 3s-orbital exchange 
hole is highly nonlocal and does not efficiently screen the 
electron, the correlation contribution goes a long way to 
make the total hole more local. 
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In addition to the obvious differences in n xc along 
each axis due to the differences in the exchange hole, 
Figs. [I] and || reveal subtle differences in the opposite- 
spin hole. The extent of the orientational anisotropy in 
the opposite-spin hole can be better visualized by plot- 
ting the hole along the x axis for the up-spin electron 
fixed on the x axis and along the z axis for the electron 
fixed on the z axis, cutting through the minima and max- 
ima of the contour plots Figs. |l](b) and||(b). These are 
represented as solid lines in Figs. ||(a) and (b). The most 
notable difference between the two cases is the height of 
the peak on the side of the atom opposite the reference 
electron, which is three times as large along the z axis (b) 
as on the x axis (a). In addition the minimum is slightly 
deeper for the z-axis case. 

Additionally, in Figs. ||(a) and (b) we show trends in 
the opposite-spin correlation hole as one gradually re- 
moves the reference electron from the atom. In addition 
to the 1.4 ag case discussed above, we place a spin-up 
electron on the x axis (a) and the z axis (b) , at a distance 
of 2.0 and 4.0 clb from atom center, and plot the on-axis 
response of the spin-down electron as long- and short- 
dashed lines respectively. These three reference radii 
correspond to placing an electron at the peak valence 
density along either axis, at the average radius from the 
atom, and at a low density point outside the atom respec- 
tively. As the electron is moved to lower densities, the 
shape of the minimum slowly gets wider and shallower, 
consistent with trends in the homogeneous electron gas. 
The position of the hole minimum stays near the atom, 
and thus increasingly more off center with respect to the 
electron. This is expected: the correlation hole, measur- 
ing the change in density in the presence of an electron 
at some reference point, can have an absolute value no 
greater than the density itself. Along the x axis, this 
trend to a shallow off-centered hole is correlated with a 
gradual increase of peak height on opposite side of atom. 
For the z-axis case the peak height remains roughly con- 
stant as the electron is removed. 

It is interesting to compare these results to recent an- 
alytic studies of the asymptotic limit of the exchange- 
correlation hole of Be and Ne ]8|,^9| . In the limit of an 
electron far removed from the atom, n xc measures the 
collapse of the remaining N — 1 electrons to an eigen- 
state of the positive ion. The correlation component of 
n xc measures the change in electron density due to the 
reduced Coulomb reduction between the remaining elec- 
trons. At intermediate distances, a few atomic radii from 
the atom center, n xc also shows a dipole polarization of 
the atom, or a reduction of electron density from the side 
of the atom nearer the reference electron to the farther 
side, more so for the more polarizable Be than for Ne. 
The collapse to the ion in the asymptotic limit results in 
a correlation hole peaked in the center with a minimum 
around the edge of the atom. This may possibly be seen 
in our data for the 4.0 as £-axis case of Fig. ||(a), where 
the correlation-hole profile develops a peak on both sides 
of the ion core. The correlation hole along both axes 
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FIG. 3. Spin-down response to a spin-up electron: (a) the 
correlation hole along the x axis for a spin-up electron fixed on 
the x axis at positions 1.4 (solid line), 2.0 (long-dashed) and 
4.0 as (dashed line) from the atom origin, (b) the correlation 
hole along the z axis for an electron fixed at corresponding 
positions on the z axis, (c) the pair correlation function for 
the 1.4 as case of (a), for energy cutoff of 20, 28 (both dotted 
lines) and 32 Ry (solid line), (d) the pair correlation function 
for the 1.4 as case of (b) for three trial wavefunctions: ips-j 
(dashed line), ifrci (dotted), tpci-j (solid line). In (a) and 
(b), the mean spin-down density is plotted for comparison 
as a dotted line. Additionally, the circles show the value for 
the "on-top" hole for each case plotted and the crosses, the 
LSD on-top hole for various points along each axis. Error 
bars show the standard deviation of the VMC expectation 
measured over several statistically independent runs. 

shows evidence of a dipole polarization of the valence 
shell when the reference electron is moved outside the 
atom, with a larger response on the "Be-like" z axis than 
along the "Ne-like" x axis. The contrast between Si and 
closed-shell atoms is that this polarization of the valence 
shell does not die out when a reference electron on the z 
axis is moved into the center of the valence shell, if any- 
thing becoming more pronounced, while the polarization 
is damped out along the x axis under the same circum- 
stances. An added complexity for Si is that the final state 
of the system after removing an up-spin electron in dif- 
ferent directions corresponds to different atomic configu- 
rations: 3s3p 2 for the z-axis case (removal of the up-spin 
3s orbital) but 3s 2 3p for the x axis. Thus significant 
differences between the correlation and exchange holes 
induced by electrons on different axes remain even for 
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electrons asymptotically far from the atom. 

We have also studied the local features of the exchange- 
correlation holes plotted, especially the opposite-spin on- 
top correlation hole, n xc (r 1 a; r, — a). This measures the 
reduction in the opposite-spin electron density at the lo- 
cation of a given electron. As the correlation hole is deep- 
est when the Coulomb repulsion is largest, the on-top 
hole also represents the lower bound for the correlation 
hole at any position in the atom | pQ ]. In previous pa- 
pers, it was observed that this feature was well described 
in several small atoms [^],p| and in Hooke's atom j3l]] by 
a local spin-density (LSDjansatz, using the on-top hole 
of the homogeneous electron gas corresponding to the 
density and magnetization at the location of the electron 
pair: 



n L J D (v,v) = n(r)( ff ^(0,n(r),C(r)) - 1). 



(11) 



Here g he9 (u,nX) is the pair correlation function of the 
homogeneous electron gas at density n, magnetization 
£ = (nt — n\)/n, and interparticle separation u. In this 
sense, it represents a link between the homogeneous elec- 
tron gas and dynamic correlations in many inhomoge- 



neous systems |32 

In Figs. ||(a) and (b) we plot the LSD on-top hole, us- 
ing homogeneous electron gas values from Ref. P3]| , at 
each point along the x and z axes (crosses). These are 
compared to on-top holes of the Slater- Jastrow wavefunc- 
tion [Eq. (||)] evaluated specifically at the location of the 
electron for the correlation holes shown (circles at 1.4, 
2.0, and 4.0 ag). These values for the Slater-Jastrow 
on-top hole are calculated directly without recourse to 
a plane-wave expansion, using VMC with a constrained 
pair random walk f34|j , and are exact within a negligible 
statistical error. The finite energy cutoffs employed for 
the rest of the data tend to lose sharp features of the 
correlation hole, in particular that of the cusp condition 
at zero interparticle separation. The resulting disagree- 
ment with the exact (for the trial wavefunction used) val- 
ues is at worst 20% and usually better. The LSD model 
agrees fairly well with the VMC data, within 10 — 20% 
for most of the points calculated here, faithfully following 
the observed trends with respect to electron position. It 
provides a fairly reasonable lower limit for the correla- 
tion hole at all points in the atom. A calculation of the 
on-top hole about the spin-down electron was also done, 
with similar agreement with theory. 

It is interesting to consider in further detail the 
anisotropy in the opposite-spin contribution to n xc at 
peak density. In the cases shown in Figs. |]|| it can be in- 
terpreted as the measure of the polarization of the down- 
spin electron density in the presence of an up-spin elec- 
tron. At first glance, one might guess that this down-spin 
polarization should roughly be invariant with respect to 
the angular orientation of the up-spin electron since one 
is measuring the response of the down-spin s-orbital to 
the angle-independent Coulomb interaction. More pre- 
cisely, one measures the correlations induced by the op- 



timized Jastrow factor which depend only on interelec- 
tron and electron-nucleus distances. The difference in the 
opposite-spin holes Figs. 0(b) and[|(b), is thus a direct 
reflection of the influence of the determinantal structure 
of the three up-spin electrons on the correlation response 
of the down-spin one. 

In the LDA, the influence of the environment on in- 
terparticle correlation is modeled by the variation with 
local density of the correlation hole, taken to be that of 
the homogeneous electron gas ||[| at the local density. 
Typically the largest absolute and relative effects in this 
model are felt at the on-top hole. At this point it seems 
to describe the observed behavior quite well, with the 
variation in density going from the maximum of the p 
orbital on the x axis to the open z axis accounting for 
the difference in the peak density on-top holes. However, 
the large discrepancy between the maxima of the peak 
density z-axis and n-axis hole, the solid lines in Fig. |jj|a 
and|^b, as well as their relatively large size, is inconsistent 
with this model. Neither does it seem easily explained by 
local gradient corrections, as the two points are located 
at a density maximum and a saddle point respectively, 
where any correction would be only second order in the 
gradient and therefore quite small. 

Within the context of a variational wavefunction calcu- 
lation, a mechanism that can introduce an orientation de- 
pendence into the correlation hole of an open shell atom 
is the three-electron "exclusion effect" or Fermi correla- 
tion ]il"| ] in configuration-interaction theory, which should 
be observable for second row atoms with less than half 
filling in the ip valence shell. For Si in the L z — pro- 
jection considered here, there is a "nondynamic" con- 
tribution to the correlation between electrons along the 
z axis due to a 3s 2 to Zp 2 z excitation that is forbidden 
for the corresponding x- or y-axis analogs by Pauli ex- 
clusion. This contribution to the correlation hole about 
an up-spin electron on the z or "unoccupied" axis, and 
the absence of the same about one on "occupied" axes 
leads to a difference in the holes observed with reference 
points on these axes. Since this particular configuration 
is within the 3p valence shell, and thus fairly close to the 
ground state energy, it is natural to expect it to play a 
significant role in the observed difference in response. 

What is striking here is that the observed orienta- 
tional dependence of n xc is obtained with a single Slater- 
determinant configuration modified by a Jastrow factor 
that depends only on interparticle distance and the dis- 
tance of particles from the core. There are no explicit 
orbital correlations in this trial wavefunction [Eq. mh]. 
The Slater-Jastrow wavefunction, coupling a Jastrow fac- 
tor with a determinant formed from an incomplete set of 
valence orbitals, in principle could have a nonzero over- 
lap with the 3pl3p x 3p y configuration allowed in the nor- 
mal exclusion-effect model, and no overlap with the ex- 
cluded excitations. However, this would be only one of 
an infinite number of configurations implicitly included in 
the Slater-Jastrow wavefunction, tied together by a "dy- 
namic" correlation factor that gives no particular varia- 
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tional weight to any one configuration. Thus it is unlikely 
that it would fully account for the observed differences 
in the response along the occupied and open axes |35| . 

On the other hand, a more general argument involv- 
ing screening can be invoked to explain the orientational 
anisotropy in our results. If the up-spin electron is placed 
on the z-axis, its exchange hole is not centered on the 
electron position, as the electron occupies a 3s orbital 
with an isotropic spatial extension throughout the atom's 
valence shell. The net effect of the electron plus its ex- 
change hole has the nature of a dipole field along the 
z axis of the atom, with an imperfectly screened nega- 
tive charge on the side of the atom nearest the electron 
and a positively charged region due to the exchange hole 
on the other side. The large peak in the down spin elec- 
tron's response to the up-spin electron can be viewed 
then as a dipole induced by the dipole field formed by 
the imperfect exchange screening of the up-spin electron. 
In the x-axis case, the exchange hole is centered and lo- 
calized on the electron; the other two electrons are in 
sp 2 orbitals blocking out the rest of the valence shell. 
Therefore the up-spin electron is efficiently screened by 
its exchange hole, and in addition, there is no "easy" 
direction for redistributing the down-spin electron den- 
sity. In such a picture the 3s 2 — > 3p 2 excitation should 
play an important but not exclusive role in the overall 
dipole response along the z axis. This effect is observ- 
able in a Slater- Jastrow wavefunction because the Jas- 
trow factor includes correlations simultaneously between 
all particle pairs. Thus, the correlation response of the 
down-spin electron to a fixed up-spin electron involves 
not only their mutual Coulomb repulsion but that of the 
other up-spin electrons as well. The inhomogeneous dis- 
tribution of these electrons about the fixed electron com- 
bined with the effect of the external potential in essence 
create a Coulomb interaction with the exchange hole. 



V. RESULTS FOR THE PAIR CORRELATION 
FUNCTION 

A. Single Configuration Wavefunction 

In a system such as an atom, the pair correlation func- 
tion g(ro, (To; r, a) is a valuable tool because, given the 
large variation in the density over the length-scale of the 
exchange-correlation hole, the shape of the hole is to a 
large extent determined by the variation in density. The 
pair correlation function can distinguish the intrinsic ef- 
fects of correlation by eliminating any features that are 
simply proportional to the density. 

In Fig. [| we plot the pair correlation function (PCF) 
for several of the cases considered previously. Specifically, 
we show g{vo, f; r, j) fixing a spin-up particle at the peak 
in the valence density, 1.4 on the x axis (a) and the 
z axis (b) and plotting the variation with respect to the 
spin-down electron in the x-z plane. These results are ob- 



tained by dividing 

fixe 

plotted in Figs. 0(b) and |(b) 
respectively by the spin-down density, with a constant 
shift of 1 due to the differing conventions in their defini- 
tions [Eqs. (0) and (|)]. The value g = 1 corresponding 
to n xc = (shown as a thick contour in Fig. 0) denotes 
the boundary between the region in which the spin-down 
density has been reduced (g < 1) and that in which it 
has been enhanced (g> 1) in the presence of the spin-up 
reference electron. As described in the next section the 
PCF can be significantly affected by statistical error at 
low density, so that the plot range is restricted to higher 
density regions where the statistical error is less than half 
a contour increment in cither direction. 

The contours of the pair correlation function in the 
vicinity of the spin-up electron on the x axis are nearly 
circular and centered on the electron, showing the spin- 
down electron repelled from it without significant direc- 
tional bias. The g — l contour is slightly oblong in shape, 
indicating that at longer distances there is a slight bias 
towards shifting electron density to the high density areas 
on the far side of the atom [x < 0). A similar plot cut- 
ting through the x-y plane (ie, cutting through both the 
p orbitals along a plane perpendicular to the one shown) 
shows similar results except that the g — l contour is more 
circular in shape. This picture is fairly consistent with 
that of a homogeneous electron gas PCF which depends 
only on the distance between particles, and indicates that 
the underlying physical processes are quite similar: high 
energy "dynamic" correlations that are quickly screened 
out over the radius of the exchange hole. 

In comparison, the PCF in Fig. 0(b), showing the re- 
sponse to the spin-up electron on the z axis shows a much 
greater departure from isotropy. The function has a no- 
ticeable dipole component with respect to the atom cen- 
ter, with a peak in g directly opposite the minimum on 
the z axis. The g—l contour has a much more shallow 
curvature than that of Fig. 0(a), so that it is no longer 
centered on the spin-up electron, but rather coincides 
roughly with the z = plane in the high density region 
of the atom. (For the 3s 2 3p x 3p y configuration of Si, the 
PCF with a reference electron on the z axis is rotationally 
invariant about the z axis so that the contours shown in 
the plot can be extended to surfaces of rotation in three 
dimensions.) In the near vicinity of the fixed electron, 
the pair correlation function is no longer a function of 
the distance between the two coordinates r and ro alone, 
but has taken on a noticeable angular dependence as well. 
Thus the slope of the well near 1.4 as on the z axis is 
steeper towards the center of the atom, shallower towards 
the outer edge of the atom. Although such a shape could 
in principle be accounted for in a region where the gra- 
dient of the density is large, using a gradient expansion 
of the homogeneous electron gas 36 , it cannot be ex- 



plained along these lines in the current situation since it 
occurs at local saddle point in the density (the peak of 
the density along the z axis) . These features in the PCF 
confirm the existence of a genuine directional anisotropy 
in the response to an electron located perpendicular to 
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(a) (b) (c) 




FIG. 4. The opposite-spin pair correlation function g(ro, T; r , I) as a function of spin-down electron position in the x-z plane 
for a spin-up particle fixed at (a) 1.4 as from the atom center on the x axis, (b) 1.4 as from atom center on the z axis, (c) 
the same as (b), with nondynamic correlations from the 3s 2 to 3p 2 double substitution included. The contour increment is 0.1 
with the contour for g— 1 represented by a thicker line; the regions between are shaded darker for lower values of g. 



the occupied p orbitals, and suggest that the explanation 
lies in a nonlocal mechanism such as the response to a 
poorly localized 3s exchange hole. 

Along this line, it is instructive to consider the 
exchange-correlation hole about the lone down-spin elec- 
tron. In this case the exchange hole is of necessity due 
to the 3s orbital regardless of the electron position, and 
one can expect significant departures from a compact, 
isotropic PCF. In the case of a down-spin electron on the 
z axis, we observe the PCF that is very close to that 
about an up-spin electron, with perhaps a slightly larger 
dipole component. In both cases, the PCF indicates the 
presence of a significant dipole response that compen- 
sates for the poorly screening 3s exchange hole. The 
PCF for a down-spin electron on the x > axis has a 
quite complicated pattern outside the on-top hole region 
- neither roughly isotropic like Fig. ||(a) or with simple 
dipole anisotropy like Fig. |^(b). There are several uncon- 
nected regions where the electron density is enhanced: 
surrounding the electron for x > and focused on the 
x axis opposite the electron for x < 0. In this case the 
exchange hole should not screen the electron efficiently, 
but there is no lowlying 3p orbital with which to con- 
struct a response, so that there is a possibility that the 
contributions of higher-order orbitals might be observ- 
able. 



B. Two Configuration Wavefunction 

Since we have observed a significant change in shape of 
the opposite-spin PCF with respect the angular orienta- 
tion of the reference electron, using a single-determinant 
Slater- Jastrow wavefunction, it is interesting to include 
explicitly the 3s 2 — > 3p 2 z excitation that is predicted in 
CI theory to play a prominent role in producing such a 
behavior in the Si atom. We consider the multideter- 



minant wavefunction ipci that consists of a 3p 2 3p x 3p y 
excited state configuration as well as the noninteracting 
ground state. This wavefunction can be treated alone or 
used as a starting point for adding further correlations 
via the Jastrow factor [Eq. (g)] ; in either case the mixing 
amplitude 771 for the excited state is variationally opti- 
mized, along with the Jastrow parameters for the latter 
case. We obtain a mixing coefficient r)i of 0.132 and total 
energy 3.7263(18) a.u. in the former case and 0.056 and 
3.8041(2) a.u. in the latter. 

The nondynamic correlation modeled by the two- 
configuration wavefunction alone, tpcii is to first order 
in rji a dipole-dipole correlation: 



5c/(ro,t;r, |) = 1 + 27ft 



z z R p (r ) Rp(r) 



r Q r R s (r ) R s {r) 



(12) 



where R p and R s are the radial 3s and 3p orbitals. The 
signature of this PCF in a contour plot with the cut 
through the x-z plane and tq fixed on the z axis is a 
series of roughly straight lines arranged antisymmetri- 
cally about the z = plane, with the g — 1 contour at 
z = 0. The shape of the function is independent of the 
position of the fixed particle, with the only change being 
in the overall amplitude. For one electron fixed on the 
x axis, on the node of the p z orbital, the PCF is to first 
order zero. 

In Fig. [|(c) we show a contour plot for the opposite- 
spin PCF of the combined CI plus Jastrow wavefunction, 
tpci-Jt f° r the same plot parameters as Fig. (|(b): with 
the up-spin electron fixed at 1.4 as on the z axis. In 
this case the occupation probability of the excited-state 
Slater- Jastrow wavefunction is very small, 0.3%, and the 
change in the total ground-state energy from the single- 
determinant Slater-Jastrow wavefunction is correspond- 
ingly small. On the other hand the explicit addition of 
the nondynamic correlation increases the correlation en- 
ergy by 1.6% and has a significant impact on the shape 
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of the hole. The PCF has a clear dipole-dipole signature, 
exaggerating the spatial anisotropy already visible in the 
one determinant case. The g = 1 contour is closer to the 
z~Q axis, and the short-range well about the fixed parti- 
cle, which for a particle near the peak of the valence shell 
density could be expected to be fairly deep and isotropic, 
has been reduced to a shallow and open-ended dip. As 
expected from the form of the nondynamic correlation, 
little change was observed in the PCF for a particle fixed 
on the x axis. 

A quantitative comparison of g for the above cases is 
also instructive and is shown in Figs. ||(c) and (d). We 
plot g(ro, T; r , i) varying r along the x axis and fixing ro 
at 1.4 as on the x axis (c) and the analogous situation for 
the z-axis (d). Error bars on the statistical measurements 
of these quantities are plotted, as well as the plane-wave 
cutoff dependence for the x-axis case. The converged 
result on the x axis (solid line) is well localized about 
the reference particle's position, with an enhancement of 
particle density on the opposite side of the atom. 

The results for the z axis, Fig. ||(d) include the PCF for 
the three different trial wavefunctions discussed above: 
the single Slater determinant plus Jastrow factor, ips~J 
(dashed line), the two-configuration wavefunction i]jci 
(dotted line) and the same multiplied by a Jastrow fac- 
tor ipci-J (solid line). In comparison to the most ac- 
curate result, ipci-j, Tpci predicts with some success 
the long-wavelength polarization response to the elec- 
tron at ro, determining how much the electron density 
is pushed from one side of the atom to the other. The 
major difference is the absence of the small dip in the 
immediate vicinity of the electron due to the cusp con- 
dition, which naturally is not obtained from the two- 
configuration wavefunction. In contrast, V's-J? which 
does not include explicit orbital-dependent correlations 
but can be optimized to obtain accurate results for short 
interparticle distances, is nearly identical to ijJci-J hi 
the vicinity of the reference electron but underestimates 
the long-wavelength polarization of the atom, obtaining 
about 70% of the PCF of ipci-J ° n the other side of the 
atom. The mixing coefficient i]i for the excited-state con- 
figuration was reduced by about 60% when the Jastrow 
factor was added to the multiconfiguration result, indi- 
cating that part of the polarization of the atom produced 
by the simple two configuration wavefunction is already 
accounted for by the Jastrow factor. In each case the 
z-axis PCF leads to a much larger peak on the far side 
of the atom than the x-axis case. 



C. Error Analysis 

There are three sources of error in our calculation: sta- 
tistical error in taking Monte Carlo estimates, the finite 
plane-wave cutoff of the data, and finally the variational 
bias due to the discrepancies of the trial from the true 
ground state wavefunction. The first two are closely con- 



nected and to some extent can be regulated; the third is 
harder to assess. 

In a typical Monte Carlo calculation, sample points 
for evaluating the density or other single-particle ex- 
pectation in a given region of space are generated 
with frequency proportional to the density itself. This 
leads to statistically precise measurements of the den- 
sity at high density and large relative errors that vary 
roughly as 1/ \/n(r) at the vanishingly low densities 
outside the atom. Calculating the VMC expectations, 
< exp(—iG -Tj) >, of a set of plane-waves periodic on 
a supercell is equivalent to taking the Fourier transform 
of a histogram distribution of statistical sample points 
on that cell. The statistical outliers from the poorly 
sampled, asymptotic low-density region show up in the 
Fourier transform as a noise background independent of 
energy. This noise can to some degree be controlled by 
imposing a finite cutoff in reciprocal space. However, 
too small a cutoff leads to spurious long wavelength os- 
cillations and is particularly poor for the short-distance 
region of the hole where the cusp condition results in a 
long range tail in the reciprocal space. We find that a 
good balance between controlling statistical and plane- 
wave error can be found by choosing a cutoff when the 
statistical average of the plane-wave component of the 
density is roughly equal to the statistical noise in its cal- 
culation. For a sampling size of around 10 5 independent 
configurations, this cutoff limit proves to be about 32 Ry 
for a resolution of 1.0 as- 

The statistical error of the Monte Carlo sampling was 
measured both for the individual plane- wave components 
of the density and conditional density, and for their real- 
space counterparts by measuring the standard deviation 
of these quantities over 10 to 20 independent runs. As 
shown in Fig. ^(c) and (d), the error bars of the PCF 
do in fact vary roughly as 1/ \/n(r) with well controlled 
errors in the peak density region, increasing to arbitrarily 
large values as one moves outside the atom. With the 
cutoff used, statistical errors in the PCF are limited to 
a value of less than 0.05 (out of a range in the PCF of 
the order of 1.0) for particles within 3.5 as of the atom 
center, which in effect provides the limits in the plots of 
the PCF in Fig. |. 

The convergence of the plane-wave expansion was 
checked by plotting the PCF as a function of cutoff en- 
ergy. A typical result is shown in Fig. ||(c) where the 
x-axis PCF for a particle fixed at 1.4 as on the axis is 
plotted, for cutoff energies of 20, 28 and 32 Ry. The exact 
on-top hole has been measured by a direct VMC calcu- 
lation and is plotted as a circle. The 20 Ry calculation 
shows clear deviations from the higher energy cutoff data, 
particularly in the region of the on-top hole, ie., in the 
immediate vicinity of the fixed electron, where the cusp 
condition contributes a high-energy tail to the exchange- 
correlation hole. The agreement between the two higher 
energy plots is well within statistical error except in the 
core region of the atom and in the low density tails. The 
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32 Ry case has not yet converged in the on-top region to 
the exact on-top value plotted as a circle, indicative of 
the slow convergence of the plane-wave expansion to the 
on-top hole cusp. 

A final source of error is from the effect of the de- 
viation of the variational trial wavefunctions [Eqs. 
and (||)] from the exact ground state. The ground-state 
energy, being variationally optimized, is typically de- 
termined with much less error than other expectations 
(though one may expect that those important in the de- 
termination of the energy, such as the density and n xc , 
should still be robust. An example of this problem is 
demonstrated in Fig. |4| in which two variationally opti- 
mized wavefunctions having relatively insignificant dif- 
ferences in total energy, give PCF's with noticeable dif- 
ferences for an electron fixed on the unoccupied z axis. 
The variational method is in general more sensitive to 
errors in the wavcfunction at high density as this region 
contributes the most to the variational energy. Thus we 
expect worse errors in the PCF when one or both coor- 
dinates are at low density. 

In order to gauge possible errors arising from varia- 
tional bias in our wavefunction, we study the change in 
the PCF and in n xc for several choices of trial wave- 
function. Specifically we used trial wavefunctions with 
four, ten and eighteen Boys and Handy basis functions, 
corresponding to keeping expansion terms of order o — 
l+m+n up to 2, 4 and 6 respectively in our Jastrow func- 
tion. The lowest order function has, in addition to two 
terms that set the cusp condition for each spin compo- 
nent, only one electron-electron and one electron- nucleus 
term. The correlation energy of each wavefunction is 
0.0640(8), 0.0812(4), and 0.0840(2) a.u. respectively, in 
comparison to 0.0877(2) a.u. using the multideterminant 
Slater- Jastrow wavefunction and 0.0883 a.u. for the CI 
calculation of Rcf. |1| . 

In Fig. H we plot the opposite-spin PCF using these 
trial wavefunctions for two cases previously considered 
in Sec. (a) an up-spin electron placed at 1.4 as on 
the x axis and (b) an up-spin electron placed at 4.0 as 
on the z axis. As in Fig. || we present a cut through 
the atom center and the reference electron, thus showing 
the minimum and maximum of the PCF. The first case 
represents a probable best-case situation, with the refer- 
ence electron placed at peak density on an axis occupied 
by a 3p orbital, so that nondynamic correlations due to 
the open-shell structure should be neglible. In contrast 
the second case presents the worst case scenario: at low 
density along the unoccupied axis. The behavior of the 
PCF in the asymptotic low-density region of the atom is 
dominated by the statistical error of our plane- wave basis 
and is not shown. 

In our best-case scenario, (a), the PCF near the refer- 
ence electron is very well obtained with even the crud- 
est model used (o = 2). On the far side of the atom 
with respect to the reference electron position, there is a 
noticeable increase in the magnitude of the peak of the 
PCF with the increase in accuracy of the wavefunction. 



r =1.4a B (100) r =4.0a B (001) 







I 

Vhf-j. 1-2 
- - - Vhf-j. 1= 4 
- - vtf HF _j, n=6 

Vci-j. n =6 


\ 

U*Xll 






(a) 

i 


(b) 


\\ 

w 

\ 

\ 


i\ 
% 

\ T I 

V X- 

X \ ' 
X 3 



-3.0 0.0 0.0 3.0 



x axis z axis 

r/a B 

FIG. 5. Convergence data for the opposite-spin PCF for 
up-spin electron fixed at (a) 1.4 as on the x axis and (b) 
4.0 as on the z axis. Cases shown are for Boys and Handy 
terms up to order o — 2,4, and 6 as described in the text, in a 
one-determinant Slater- Jastrow wavefunction, and o — 6 for a 
two-determinant Slater- Jastrow wavefunction. Additionally, 
results for a two-determinant wavefunction without Jastrow 
factor are given in (b) as a dot-dashed line. 

As expected, the addition of nondynamic correlations in 
form of a two-determinant Slater- Jastrow wavefunction 
(solid line) has no noticeable effect. An investigation of 
the corresponding contour plots, cutting through the x-z 
plane as shown in Fig. ||, shows that the depth and spa- 
tial extent of the PCF is obtained with the lowest order 
wavefunction and varies only slightly with increase in ba- 
sis set. The major difference is the gradual adjustment 
of the 5 = 1 contour from a more isotropic shape to the 
elliptical one shown in Fig. ||(a), which accounts for the 
gradual increase in the peak on the far side of the x axis 
as the g = 1 contour moves towards the atom center. 

The worst case scenario (b) shows a far greater degree 
of disagreement between wavefunctions. The three one- 
determinant Slater- Jastrow wavefunctions agree fairly 
well with each other, with an increasing percentage of 
the electron density removed from the near side of the 
atom and placed on the far side. The peak height on 
the far side is roughly the same as for the 1.4 as, high- 
density case shown in Fig. |^(d) . The introduction of the 
nondynamic 3s 2 — > 3p 2 excitation into the Slater- Jastrow 
wavefunction (solid line) leads to a dramatic change in 
the shape of the PCF. The peak of the PCF from the 
base value of g = 1 representing no change due to corre- 
lations is increased by a factor of three, and only at the 
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on-top hole (not shown) is there no significant change 
in the PCF. In comparison, the PCF of the optimized 
two-determinant wavefunction ipci discussed in Sec. [v]is 
shown as a dot-dashed line. The two-determinant Slater- 
Jastrow wavefunction is an interpolation between the two 
limiting cases, favoring tpci at most locations. 

An explanation of the "gigantic" features in the non- 
dynamic part of the PCF in Fig. |5|(b) comes from the 
observation that the ratio zR p (r)/R s (r) between the 3p 
and 3s orbital, which determines the shape of the nondy- 
namic PCF to first order [Eq.|l§, increases exponentially 
at large distances along the z axis. As a result, in the 
asymptotic region (that is, either the response to a refer- 
ence electron outside the atom or the asymptotic tails of 
the response to a reference electron at high density), the 
nondynamic correlation introduced with the 3s 2 — > 'ip 2 z 
substitution is no longer a small perturbation even with 
a small mixing parameter. The Jastrow factor basis func- 
tions however are chosen to tend to a constant at either 
large electron-electron or electron-nucleus distances. It 
is quite likely that none of the cases studied accurately 
represents the true asymptotic behavior of the correlation 
function. They should rather be considered to provide a 
qualitative idea of the PCF as well as a sense of the range 
of behavior it should reasonably lie within. 

It is interesting to note that the agreement between 
the various PCF's is far closer for the high density case 
on the z axis plotted in Fig. ||(d). This indicates the 
greater weight of the exchange-correlation hole at high 
density in determining the total exchange-correlation en- 
ergy, and the corresponding robustness of its determi- 
nation even with qualitatively different wavefunctions. 
Given the large degree of variation in the asymptotic 
limit, it would be interesting to study the effect of in- 
cluding a larger number of configurations in a multidc- 
tcrminantal Slater- Jastrow wavefunction or include the 
effects of backflow or multielectron coordinates Jlij into 
the evaluation of orbitals. At higher densities we ex- 
pect the effect of such improvements in the wavefunction 
would be to add further refinements in the shape of the 
PCF along the z axis, but with much less change in basic 
features such as its range and magnitude. 

VI. CORRELATION EFFECTS ON THE SPIN 
DENSITY 

In addition to the exchange-correlation hole, the spin 
components of the single-particle density change with 
the inclusion of Coulomb correlations. In principle, the 
Kohn-Sham equations used to derive the density in DFT 
should give the exact ground-state radial density, even for 
a degenerate ground-state |T^| ; however there is no such 
principle for the spin components of the density. Thus 
a change observed in the spin components of the density 
that does not alter the total radial density can be con- 
sidered an intrinsic feature of Coulomb correlation and 



not simply due to the inaccuracy of the LDA density. As 
with n xc , a prominent feature of the change in the spin- 
dependent density, An(r,er), is anisotropy with respect 
to the fully occupied x and open z axes. This anisotropy 
in An(r, a) may help to shed light on the mechanisms 
underlying the anisotropy in n xc . 
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FIG. 6. The change in the spin-dependent density 
ni(r, a) — no(r, er) between the fully interacting and noninter- 
acting systems. The left side of the plot shows a cut through 
the x axis; the right side, through the z axis. The solid line 
is the spin-down density change using the one-determinant 
Slater- Jastrow wavefunction ips-j', the long-dashed is the 
spin-up density change for the same wavefunction. Dot- 
ted and dashed lines give the corresponding quantities for 
the two-determinant Slater- Jastrow function ipci-j- The 
dot-dashed line is for the two-determinant function ipci- The 
plot units are 10~ 3 a.u. 

In Fig. | the spin-down density for the single determi- 
nant Slater-Jastrow wavefunction is plotted as a dotted 
line along the x axis (a) and z axis (b) . Note that in the 
nonintcracting wavefunction this density is that of the 
3s orbital and thus the same along each axis. However, 
in the interacting case, anisotropy in the correlation re- 
sponse reduces the density along the x axis by roughly 
6% and increases it by the same amount along the z- 
axis. This change is shown in Fig. ^ (solid line) along 
with the change in the spin- up density (long-dashed line) , 
with the right side of the plot showing a cut through the 
x axis and the left side a cut through the z-axis. The 
corresponding density changes using the two-determinant 
Slater-Jastrow wavefunction (ipci~j) are plotted as dot- 
ted and dashed lines. In both cases, the change in the 
spin-up density has a qualitative trend opposite to that 
of the spin-down density, showing a large enhancement 
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in density near the peak of the 3p x orbital and a drop 
in density along the z axis. Thus, the changes in the 
two spin components cancel out partially in the total 
density, but add to the total spin density, defined as 
m(r) = n(r,1) — n(r, 1). In addition there are observ- 
able effects of a small reduction in the radius of the atom 
that occurs with the addition of the Jastrow correlation 
term, particularly a reduction in up-spin density at a 
large distance from the atom and an enhancement of the 
down-spin density in the core region. A convergence test 
along the lines of that for n xc in the previous section was 
carried out for the spin components of the density. The 
qualitative trends are repeatable with the less accurate 
Jastrow factors but with considerably slower quantitative 
convergence than for n xc . 

An anisotropic shift in spin density can be gener- 
ated from the nondynamic Fermi correlation discussed in 
the previous sections: the mixing of the noninteracting 
ground-state with a 3s 2 — > 3p 2 -substitution excited state 
and the lack thereof due to Pauli exclusion for the 3p x 
and 3p y analogs. Using the two-determinant wavefunc- 
tion ipci, the resulting change in each spin component of 
the density is 

An(r,T) = An(r,|) =r? 2 (|^(r)| 2 - |fe(r)| 2 ) (13) 

where ip3 S and ip3 Pz are the wavefunctions of the 3s and 
3p z orbitals and r/i the excited state probability ampli- 
tude. This function, shown as a dash-dotted line in Fig||, 
corresponds to a density enhancement along the z axis 
and reduction in the x-y plane for both spin-components 
for a net zero change in the spin density m(r), in marked 
contrast to the mutually opposing changes of the two 
spin-components of the Monte Carlo data. Given the op- 
timal value of rjf = 0.017 for the occupation number of 
the excited-state configuration, one finds that the magni- 
tude of the change in either spin component of the den- 
sity is much smaller than observed, and qualitatively in 
the wrong direction for the up-spin case. When the non- 
dynamic s 2 — » p 2 z excitation is included into the Slater- 
Jastrow wavefunction (dotted and dashed lines), both the 
up-spin and down-spin densities are enhanced slightly 
along the z axis, but the qualitative picture remains 
unchanged - as one might expect from the very small 
value of r/f = 0.003 that was optimal for this case. For 
these reasons, it is unlikely that the observed anisotropic 
change in the spin components of the density can be ex- 
plained by this type of mechanism. 

On the other hand, a simple explanation of these re- 
sults lies in that the observed change in the spin density 
reduces the spatial overlap between the two spin compo- 
nents. As the correlation energy is predominantly deter- 
mined by the spatial correlation between opposite-spin 
electrons, such a change in spin density can reduce the 
correlation energy in a way that reduces the total energy 
if the total density remains unchanged. This correlation 
effect is similar to that of the unrestricted Hartree-Fock 
method pTj in which the energy of a system like H2 can 



be lowered from its Hartree-Fock value by breaking a 
symmetry of the ground state to induce the spatial sep- 
aration of opposite spin-components of the density. In 
the Si atom, since the ground state is degenerate and 
lacks the symmetry of the Hamiltonian, it already has 
a nonzero spin density; multiplication by a Jastrow fac- 
tor does not break the symmetry of the ground state 
in a substantial way [3q| . The change in spin-density re- 
flects the effect of the Coulomb interaction which induces 
the spatial separation between opposite spin electrons, in 
a system in which the different spin components are to 
some degree already spatially separated in the noninter- 
acting state. In contrast to a fillcd-shell atom or other 
spin-unpolarized system, this spatial separation of oppo- 
site spins appears not only as a correlation hole but as a 
change in the mean density as well. Consistent with this 
picture, the direction of the change in the L z = projec- 
tion of the Si atom enhances the absolute difference in 
the spin density n(r, f) — n(r, |)| that already exists in 
the Hartree-Fock wavefunction, where up-spin electrons 
occupies p x and p y orbitals while the lone down-spin elec- 
tron does not. 



VII. SUMMARY AND CONCLUSION 

We have calculated the exchange-correlation hole of 
the valence shell of the ground state of the Si atom as 
a function of spin decomposition, using the Variational 
Monte Carlo method. This relatively simple four electron 
system, restricted to a single valence shell, nevertheless 
shows a rich variety of phenomena in exchange and cor- 
relation not present in closed-shell atoms. 

The incomplete filling of the valence shell of the open- 
shell atom leads to dramatic anisotropy and nonlocality 
in the exchange hole, even at peak densities in the valence 
shell. This is accompanied by a significant compensat- 
ing anisotropy in the correlation hole making the total 
exchange-correlation hole more (but not completely) lo- 
cal and isotropic. 

Our paper has focused mostly on the exchange- 
correlation hole about a majority (up) spin electron. In 
this case, as one goes from a reference point at peak den- 
sity along an "occupied" axis, along which one of the 
occupied 3p orbitals is oriented, to one on the "open" 
z axis perpendicular to the 3p orbitals, the exchange hole 
changes from an efficiently screening 3sp 2 -like to a poorly 
screening 3s character. As a result, we observe in the re- 
sponse of the minority (down) spin density a significant 
"dipole" shift or shift of density from one side of the 
atom to the other, that occurs when the up-spin electron 
is placed in the peak density position on the open axis 
and not when it is placed in the center of a p orbital. 
This difference shows up in the pair correlation function 
as a difference in shape, with that on the "occupied" 
axis being a modest distortion of the isotropic shape of 
a dynamic correlation, and that on the z axis showing a 
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marked dipole component along the z axis. The dipole 
polarization observed is also notable in that it occurs not 
only for a reference electron outside the atom but at the 
peak along the z axis, where the local gradient is zero; 
it is therefore a truly nonlocal feature not amenable to 
modeling by an expansion in the local density gradient. 

The explicit inclusion of nondynamic correlations into 
our wavefunction enhances the difference between open 
and occupied-axis response, particularly when the refer- 
ence electron is moved outside the atom. This is due 
to the exclusion effect in which the 3s 2 — > 3p 2 substitu- 
tion is allowed along the open axis and excluded along 
the occupied ones. Nevertheless, the existence of such a 
difference in the single-determinant Slater-Jastrow wave- 
function, given the structure of the Jastrow factor, seems 
better explained in terms of the screening of the exchange 
hole, in which the Coulomb interaction of the down spin 
electron with the up-spin 3s exchange hole causes the 
dipole response along the open axis. In this case, it is the 
lowlying 3p orbital that best compensates for the poor 
screening of the exchange hole. It should thus play an 
important if not exclusive role in inducing the anisotropy 
of the opposite spin correlation hole. 

A second major effect of the open-shell structure of 
the atom upon its response to the Coulomb interaction 
is the change of the spin components of the density with 
the addition of correlations. The down-spin density is 
pushed inwards and onto the unoccupied axis and the 
up-spin density pushed off the z axis and onto the peak 
of the p orbitals in the x-y plane, resulting in a reduc- 
tion of the spatial overlap between the two spin compo- 
nents in a way which leaves the radial density largely 
unchanged. The standard exclusion effect resulting from 
the 3s 2 — > 3p 2 substitution produces a net change in den- 
sity that is both quantitatively too small and qualita- 
tively incorrect. However, both the anisotropic features 
of the exchange-correlation hole and the changes in spin 
density essentially stem from the same effect: the ten- 
dency of the Coulomb interaction to induce the spatial 
separation of opposite spin electrons in the context of the 
spatially anisotropic and spin-polarized structure of the 
degenerate Si ground state. 

We have only briefly discussed the n xc for a down-spin 
electron, where in addition to the screening of the ex- 
change hole (always 3s-like) , the anisotropy of the deter- 
minantal structure of the up-spin electrons plays an im- 
portant role. The correlation response to the down-spin 
electron on the occupied axis combines a compensation 
for poor screening by the exchange hole and the absence 
of a lowlying 3p-orbital component from the correlation 
hole on account of Pauli exclusion. We have observed 
in this situation subtle structural properties in the PCF 
that merit further investigation. 

In contrast to longer ranged features of the correlation 
hole we find that its on-top value - the reduction of elec- 
tron density in the immediate vicinity of an electron - 
is reproduced by a local density ansatz within 10-20%, 
over a fairly wide range in density and magnetization, 



and is otherwise insensitive to the structure of the atom. 
Given the energetically reasonable shape of the LDA n xc 
and that it satisfies the short-range cusp condition and 
global particle-sum rule of the true n xc in addition to 
the approximate fit to the Si on-top hole, it is reason- 
able to expect that it will provide a good approximation 
for the angle- or system-averaged n xc , that is, after av- 
eraging out many of the subtle angle or position depen- 
dent features studied here ||. Nevertheless, the complex 
spin-dependent phenomena observed in this paper point 
to the inherent difficulty of systematically improving on 
local or semilocal density functional theories in systems 
with nontrivial valence-shell structure, such as open shell 
atoms or multiply-bonded molecules. 

Our results provide support for several recent hybrid 
approaches to DFT. The combination of a short ranged 
on-top hole that is fairly well modeled by local density 
functional theory and longer-ranged exchange and exclu- 
sion effects that are not indicates the usefulness of den- 
sity functional schemes which include the on-top hole as 
a basic component jj2j or involve the hybridization of 
short-ranged local or semilocal density functionals with a 
more accurate treatment of longer ranged correlations us- 
ing RPA |36| or CI |Q . The tendency of the correlation 
hole to cancel out the anisotropy in the exchange hole, 
and increase the sensitivity of the exchange-correlation 
hole to electron position, along with the quality of the 
on-top hole in the LSD approximation, lends support to 
recent hybrid approaches p2[ which mix the exact ex- 
change hole with an local DFT approximation of the 
full-coupling constant exchange-correlation hole. These 
methods rely on the assumption that the full coupling 
constant limit of the integral in Eq. ([I]) is more likely to 
be amenable to approximation by the isotropic, localized 
n xc of the homogeneous electron gas than the noninter- 
acting limit dominated by exchange. In the case of the 
open-shell atom exclusion effects quite effectively can- 
cel out the gross anisotropy in the exchange hole caused 
by the open-shell structure, at least within the region of 
peak valence density. 

On the other hand, the changes in spin-density that 
we observe indicate that a standard Hartree-Fock or re- 
stricted CI basis set may be a less desirable starting point 
for implementing hybrid DFT methods in spin-polarized 
systems than, for example, a generalized Hartree-Fock 
approach that would allow for anisotropic distortions in 
the spin density due to Coulomb correlations. Also, it 
seems possible that more could be done to obtain ac- 
curate correlation energies within the local spin-density 
approximation with the incorporation of projection spe- 
cific information. 
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